* Evaluate recovered beliefs from GSU gender and competitiveness beliefs */

* log file
capture log close _all
log using "Evaluate Recovered Beliefs on GSU Gender and Competitiveness.log", replace name(recover)

* configure Stata
capture: version 18
set processors 4
set more off
set scheme s1color
set seed 987654321
timer clear 1
timer on 1

* tasks
global doIND "y"
global doPOOL "n"
global doREPORT "n"

* number of subjects
local Nsub = 88

* number of questions
local Nq = 2

* name of file stem
local belief "gc"

* risk model (eut or edu)
local risk_model "rdu"

* graphics font
graph set window fontface "Candara"

* create directories needed
capture: mkdir figures

* tell us what version ran
about

* evaluate over subjects and questions
if "$doIND" == "y" {

qui {
noi: di "Evaluating beliefs of `Nsub' subjects..."

* figure counter
local f = 2

* main calculations
*forvalues s=1/`Nsub' {
forvalues s=2/3 {
forvalues q=1/`Nq' {

local bin = 4

* check for file
capture: confirm file recovered/`belief'_s`s'_q`q'.dta
di _rc

* file exists
if _rc == 0 {

* get the recovered beliefs
use recovered/`belief'_s`s'_q`q', clear
summ

* risk preferences
summ r
local r = r(mean)
local r_   = string(`r', "%3.2f")
if "`risk_model'" == "rdu" {
	summ eta
	local eta = r(mean)
	local eta_ = string(`eta', "%3.2f")
	summ phi
	local phi = r(mean)
	local phi_ = string(`phi', "%3.2f")
}

* label for mcmc sample
rename mcmc mcmc_size
generate long mcmc = _n
label variable mcmc "MCMC sample number"

* check beliefs sum to 1
generate check_belief = 0
forvalues x=1/`bin' {
	replace check_belief = check_belief + b_`x'
}
summ check_belief
rename check_belief normalize_b
label variable normalize_b "Normalization for recovered beliefs to sum to 1"

* get posterior predictive mean and median for each bin
generate sum_p_mean = 0
generate sum_p_median = 0
forvalues x=1/`bin' {
	summ b_`x', detail
	local m = r(mean)
	generate p_mean_`x' = `m'
	label variable p_mean_`x' "Posterior predictive mean for bin `x'"
	local m = r(p50)
	generate p_median_`x' = `m'
	label variable p_median_`x' "Posterior predictive median for bin `x'"
	replace sum_p_mean = sum_p_mean + p_mean_`x'
	replace sum_p_median = sum_p_median + p_median_`x'

	label variable r_`x' "Report for bin `x'"
	label variable b_`x' "Recovered belief for bin `x'"
}
* normalize posterior mean and median
forvalues x=1/`bin' {
	replace p_mean_`x' = p_mean_`x'/sum_p_mean
	replace p_median_`x' = p_median_`x'/sum_p_median
}
summ p_mean_*
summ p_median_*

rename sum_p_mean normalize_mean
label variable normalize_mean "Normalization for recovered posterior predictive mean to sum to 1"
rename sum_p_median normalize_median
label variable normalize_median "Normalization for recovered posterior predictive median to sum to 1"

* illustrate
des mcmc r_* b_* normalize_b normalize_mean p_mean_* normalize_median p_median_*
list mcmc r_* r eta phi in 1/5, noobs
list mcmc b_* normalize_b in 1/5, noobs
list mcmc normalize_mean p_mean_* in 1/5, noobs
list mcmc normalize_median p_median_* in 1/5, noobs

* save, reduce to one observation and save the summary
save recovered/`belief'_s`s'_q`q'_save, replace
keep if _n == 1
save recovered/`belief'_s`s'_q`q'_summary, replace

* display
list, noobs
reshape long r_ b_ p_mean_ p_median_ ml_ true_ post_, i(mcmc) j(bbin)
list, noobs
rename r_ report
rename b_ belief
rename p_mean_ p_mean
rename p_median_ p_median
rename post_ post

* comparison of reports and beliefs
local angle = 0
if "`risk_model'" == "eut" {
	if "`q'" == "1" {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the Piece Rate Rank Question", size(medlarge)) subtitle("Posterior mean EUT risk preference is r = `r_'", size(small)) saving(figures/`belief'_s`s'_q`q', replace)
	}
	if "`q'" == "2" {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the Tournament Rank Question" , size(medlarge)) subtitle("Posterior mean EUT risk preference is r = `r_'", size(small)) saving(figures/`belief'_s`s'_q`q'_, replace)
	}
}
if "`risk_model'" == "rdu" {
	if "`q'" == "1" {
		if `s' == 2 | `s' == 3 {
			local f = `f'+1
			graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Figure C`f': Beliefs of Subject #`s' for the Piece Rate Rank Question", size(medlarge) span) subtitle("Posterior mean RDU risk preferences are r = `r_', {&eta} = `eta_' and {&phi} = `phi_'", size(small) span) saving(raven_C`f'.gph, replace)
		}
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the Piece Rate Rank Question", size(medlarge)) subtitle("Posterior mean RDU risk preferences are r = `r_', {&eta} = `eta_' and {&phi} = `phi_'", size(small)) saving(figures/`belief'_s`s'_q`q', replace)
	}
	if "`q'" == "2" {
		if `s' == 2 | `s' == 3 {
			local f = `f'+1
			graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Figure C`f': Beliefs of Subject #`s' for the Tournament Rank Question", size(medlarge) span) subtitle("Posterior mean RDU risk preferences are r = `r_', {&eta} = `eta_' and {&phi} = `phi_'", size(small) span) saving(raven_C`f'.gph, replace)
		}
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "Best" 2 "Second Best" 3 "Third Best" 4 "Fourth Best") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the Tournament Rank Question" , size(medlarge)) subtitle("Posterior mean RDU risk preferences are r = `r_', {&eta} = `eta_' and {&phi} = `phi_'", size(small)) saving(figures/`belief'_s`s'_q`q'_, replace)
	}

}
graph export figures/`belief'_s`s'_q`q'.png, replace

* close loop over skip for missing questions
}

* close loop over questions
}

* close loop over subjects
}

* end of qui
}

* end of $doIND
}

* collate files
if "$doPOOL" == "y" {

* qui block
noi {

* keep track of files -- summary files
local f = 0
forvalues s=1/`Nsub' {

	forvalues q=1/`Nq' {

		* check for file
		capture: confirm file recovered/`belief'_s`s'_q`q'_summary.dta
		di _rc

		* file exists
		if _rc == 0 {

				* get the recovered and processed beliefs
				use recovered/`belief'_s`s'_q`q'_summary, clear

				* accumulate
				local f = `f'+1
				save tmp`f', replace

		}

	}
}

* append
use tmp1, clear
forvalues x=2/`f' {
	append using tmp`x'
	erase tmp`x'.dta
}
erase tmp1.dta
save `belief'_`risk_model'_summary, replace
noi: di "* saved the pooled summary file..."

* keep track of files -- save files
local f = 0
forvalues s=1/`Nsub' {

		forvalues q=1/`Nq' {

		* check for file
		capture: confirm file recovered/`belief'_s`s'_q`q'_save.dta
		di _rc

		* file exists
		if _rc == 0 {

				* get the recovered and processed beliefs
				use recovered/`belief'_s`s'_q`q'_save, clear

				* accumulate
				local f = `f'+1
				save tmp`f', replace

		}

	}
}

* append
use tmp1, clear
forvalues x=2/`f' {
	append using tmp`x'
	erase tmp`x'.dta
}
erase tmp1.dta
save `belief'_`risk_model'_save, replace
noi: di "* saved the pooled save file with all MCMC samples..."

* end of qui
}

* review
tab sid
tab qid

* end of $doPOOL
}


* PDF document
if "$doREPORT" == "y" {

* quietly
noi {

* start a PDF
	noi: di "* Generating PDF report..."

	putpdf begin, font("Candara",12) margin(left,1) margin(right,1)

	putpdf paragraph, font("",14) halign(center)
	putpdf text ("Recovered Beliefs about Gender and Competitiveness"), bold linebreak(2)

	local d = c(current_date)
	local t = c(current_time)
	putpdf paragraph, font("",12) halign(center)
	putpdf text ("`d'   `t'"), linebreak(2)

	putpdf paragraph, font("",12) halign(left)

	putpdf text ("        There have been `Nsub' subjects, each reporting their beliefs about their rank in terms of their Piece Rate score and their Tournament score. ")
	putpdf text ("These were tasks 5 and 6, respectively, of the GSU replication and extension of the Niederle and Vesterlund study of gender and competitiveness. ")
	putpdf text (" "), linebreak(2)

	putpdf text ("        The belief interface we used employs a Quadratic Scoring Rule (QSR) ")
	putpdf text ("with parameters chosen to allow a perfect allocation of 100 tokens to the correct interval to be rewarded with $30 if that questions was chosen for payment. ")
	putpdf text ("The QSR reports were defined over 4 bin intervals, whether the rank was First, Second, Third or Fourth. ")
	if "`risk_model'" == "eut" {
		putpdf text ("The risk preferences were estimated with a Bayesian Hierarchical Model (BHM), assuming an EUT model of risk preferences. ")
	}
	if "`risk_model'" == "rdu" {
		putpdf text ("The risk preferences were estimated with a Bayesian Hierarchical Model (BHM), assuming a RDU model of risk preferences with Prelec probability weighting. ")
	}
	putpdf text ("Beliefs were recovered from these reports using the individual risk estimates from the BHM for each subject, and assuming that they responded rationally when reporting beliefs conditional on those risk preferences. ")
	putpdf text (" "), linebreak(2)

	* list the last page of text
	local page = 2

	* individual beliefs
	forvalues s = 1/`Nsub' {

		putpdf pagebreak
		putpdf paragraph, halign(right) font("Candara",10)
		local page = `page' + 1
		putpdf text ("Page `page'"), linebreak(1)
		putpdf paragraph, halign(center)
		putpdf image figures/`belief'_s`s'_q1.png, linebreak width(6)

		* check for file -- one subject did not do Task 6
		capture: confirm file figures/`belief'_s`s'_q2.png
		* file exists
		if _rc == 0 {
			putpdf image figures/`belief'_s`s'_q2.png, linebreak width(6)
		}

	}

* save the PDF report
	if "`risk_model'" == "eut" {
		putpdf save "Recovered Beliefs about Gender and Competitiveness using EUT.pdf", replace
	}
	if "`risk_model'" == "rdu" {
		putpdf save "Recovered Beliefs about Gender and Competitiveness using RDU.pdf", replace
	}

* end of quietly
}

* end of $doREPORT
}

* time taken overall
timer off 1
timer list
local secs = r(t1)
local mins = `secs'/60
local hrs = `mins'/60
local secs_ = string(`secs', "%10.0f")
local mins_ = string(`mins', "%4.1f")
local hrs_ = string(`hrs', "%4.2f")
di "Complete calculations took `secs_' seconds, `mins_' minutes, or `hrs_' hours."

log close recover
